	PRO TVS_JZ,dir,ftime,Nproc,jj,OP=OP,PS=PS
; Interactive subroutine
; You have to set the X-point with cursor function
; (Halftones of the electric current density component along the z-direction, j(x,y))
; If OP=1 then you have to set the O-point too

leer_input,dir,N,nu,eta,eps,epse,dt,nstep,nstg,nstdat,nstk,nstp
loadct,1 ; color palette
time = '000'+string(ftime,format='(i1)')
if ftime gt 9 then time = '00'+string(ftime,format='(i2)')

x=2.*!pi*dindgen(N)/N & y=x
t = 0.0d0
apX=dblarr(ftime+1) & apO=apX  & recflux=apX & tiempo=apX 

; For each frame (time)
for i=0,ftime do begin
time = '000'+string(i,format='(i1)')
if i gt 9 then time = '00'+string(i,format='(i2)')
; Compute potencial a
leer,dir,time,Nproc,a,f,b,u,jj,ww,t
tiempo(i) = t/eps
window,0,xs=N,ys=N,title='j(x,y)'
tvscl,jj
cursor,xpx,ypx,/device,/wait
print,'Frame: ',i
apX(i) = a(xpx,ypx) ; guarda el valor de a en el punto x
apO(i) = a(N/2-1,3*N/4-1)
endfor

if keyword_set(OP) then begin
print,'======================================================================='
print,'PUNTO-O'
print,'======================================================================='
t = 0.0d0
for i=0,ftime do begin
time = '000'+string(i,format='(i1)')
if i gt 9 then time = '00'+string(i,format='(i2)')
; Compute potencial a
leer,dir,time,Nproc,a,f,b,u,jj,ww,t
tiempo(i) = t/eps
window,0,xs=N,ys=N,title='j(x,y)'
tvscl,jj
cursor,xpo,ypo,/device,/wait
print,'Frame: ',i
apO(i) = a(xpo,ypo) ; guarda el valor de a en el punto O
endfor
endif

apX = apX(0) - apX & apO = apO(0) - apO 
recflux = (apX - apO)/eps
; Plot
if not keyword_set(ps) then window,0,xs=600,ys=600
if keyword_set(ps) then begin
set_plot,'PS'
device,filename=dir+'recflux_'+time+'.eps',bits_per_pixel=8,/encapsulated,xs=15,ys=12
endif
loadct,0 ; palette black-white
plot,tiempo,recflux,xtitle='Time',ytit='Reconnected Flux',xr=[0,max(tiempo)],yr=[0,1.1*max(recflux)],xsty=1,ysty=1,charsize=2
oplot,tiempo,recflux,psym=4
if keyword_set(ps) then begin
device,/close
set_plot,'X'
endif
end
